--- title: Graphics keywords: fastai sidebar: home_sidebar summary: "This module encapsulates the functions used for displaying an optimisation function" description: "This module encapsulates the functions used for displaying an optimisation function" nb_path: "01_graphics.ipynb" ---
The idea of zooming is actualy closely related to the min and max boundaries set on the plots:
Since zoom is usually a proportion of the current view, there's a bit of math to do in order to get the computations correct.
This section is highly specific and beside the main point but I like to keep it just because it's time spent that I don't want to be lost to history. If you're reading this, you might probbably want to skip it since it isn't clever nor informative..
So the main idea of the zooming that I'm going to implement is described by the image bellow:
mean (center) between themmin and max closer to the mean by the same amunt (percentage), f called a zooming factor.assert zoom((-5, 5), (-5, 5), 0.1) == ((-4.5, 4.5), (-4.5, 4.5))
Funny enough, these can be vectorized by the numpy transformations bellow
d = np.array([(-5, 5), (-5, 5)])
zoom_factor = 0.1
means = d.mean(axis=-1) # computing the means
distances = np.abs(d - means)
change = distances * zoom_factor # compute the change need
change_with_direction = change * [1, -1] # add signs for the direction of changes (mins should increase, maxes should decrese in value)
zoomed_d = d + change_with_direction
d + np.abs(d - d.mean(axis=-1)) * [1, -1] * zoom_factor # single line transformation
Now, we see that the above computation does, an abs where we eliminate the sign, and right after that we add it back by multipligin with [-1, 1]. If we think in terms of moving from the mean left and right we can simplify the formula a bit, as follows:
d.mean(axis=-1) - (d.mean(axis=-1) - d) * (1 - zoom_factor)
If we analitically decompose the (1 - 0.1) term and simplify the result, we remain with the bellow formula:
d + (d.mean(axis=-1) - d) * zoom_factor
Which is almost identical with the one we've started from but, as we've observed, does not trim the signs and adds them back.
from optimisations.functions import himmelblau
fig = plt.figure(figsize=(13, 10))
ax = plt.axes(projection='3d')
plot_function_3d(himmelblau(), ax=ax, azimuth=20, angle=225)
Let's make it interactive so we can play with it a bit.
Note: This code will probably not work in the blog post
import ipywidgets as widgets
from ipywidgets import interact, interact_manual
@interact
def plot_interactive(azimuth=45, rotation=45, zoom_factor=(-1, 1, 0.1), show_projections={True, False}):
return plot_function_3d(function=himmelblau(), azimuth=azimuth, angle=rotation, zoom_factor=zoom_factor, show_projections=show_projections)
Note that we have a 3D plot what can handle rotations, we need to allow this capability to the 2D plot as well, since we wish the two charts to move in sync.
What we need to do is use linear algebra and rotate the initial meshgrid of points. We accomplish this by simply multiplying it with a rotation matrix.
Adding rotation to the second plot
angle = 40
xx, yy, zz = contour(himmelblau())
radians = angle * np.pi/180
# counter-clockwise rotation matrix
# https://stackoverflow.com/questions/29708840/rotate-meshgrid-with-numpy
rotation_matrix = np.array([
[np.cos(radians), -np.sin(radians)],
[np.sin(radians), np.cos(radians)]
])
xx, yy = np.einsum('ji, mni -> jmn', rotation_matrix, np.dstack([xx, yy]))
plt.contourf(xx, yy, zz, levels=1000, cmap='Spectral', norm=colors.Normalize(vmin=zz.min(), vmax=zz.max()), alpha=0.3)
We've successfully rotated the contour, but we see that the axis were left unchanges. So we didn't rotate the plot, but merely it's contents.
At this point, I guess there are three options:
PIL, imagemagick), export the 2D plot as an image, rotate it and then display itThe first option looks like the biggest hack of all since it involves adding at least 2 new dependencies for this sole purpose (plus the additional computation).
The last one is easy to see without further experiments, so we only need to see how option 2 looks like.
Drawing the 2D plot on a 3D Axis and viewed from above. This makes it possible to dispay the rotated axis as well but also makes the plot smaller.
fig = plt.figure()
ax_ = fig.add_subplot(1, 2, 2, projection='3d')
angle = 125
## plot_function_2d
xx, yy = np.meshgrid(
np.linspace(-5, 5, num=100),
np.linspace(-5, 5, num=100)
)
zz = himmelblau()(xx, yy)
ax_.contour(xx, yy, zz, levels=200, offset=0, cmap='Spectral', norm=colors.Normalize(vmin=zz.min(), vmax=zz.max()), alpha=0.5, zorder=1)
ax_.view_init(90, angle)
plt.tight_layout()
I think I will stick with the 2D approach as the contour is bigger.
The rotate function bellow looks way more complex than we've sketched it out above, and this is because it is now handling two use cases that it needs to disambiguate:
The np.einsum part is the core of the function, and it's operation is more complex, leaving its explanation for a future post.
Another interesting bit of it is the way the code bellow works, namely how the parameters are returned.
xx, yy = np.einsum('ij, mnj -> imn', rotation_matrix, np.dstack([xx, yy]))
Normally, any np. prefixed operation returns a np.ndarray but in this case we see that the enisum function is able to destructure the return into two separate parameters xx and yy. This happens with the help of the np.dstack function.
Let's see what it does, bellow:
print(np.dstack([xx, yy]).shape)
np.dstack([xx, yy])[:3, :3, :]
So dstack just made a stack of it's arguments, on a new, rightmost axis (the last 2 value from the shape). This shape, coupled with the specified einsum operation, where both i and j are equal to 2 gives a result with the shape of (2, 100, 100) semantically similar to a tuple (xx.shape(100, 100), yy.shape(100, 100)). This shape enables the python language to destructure this return into 2 independent values, the ones we actually wanted.
The gist of single points rotation is a matrix multiplication but since the meshgrid is implemented in einsum notation as is rather elegant, at least because of the arguments decomposition trick, I'll try experimenting a bit to replicate the matrix multiplication it with einsum as well.
angle = 225
# apply rotation
angle = (angle) % 360
radians = angle * np.pi/180
# clockwise rotation matrix
rotation_matrix = np.array([
[np.cos(radians), -np.sin(radians)],
[np.sin(radians), np.cos(radians)]
])
points = np.array([
[1, 1],
[-2, 2],
[3, -3],
[4, 4]
])
function = himmelblau()
coords = function.coord(points)
coords[:, [0, 1]] @ rotation_matrix
Ok, this is the (correct) result we're getting using the plain matrix multiplication. Our goal is the einsum operation that computes these, but reshapes them as well in (2, 4) so we can later on decompose them as (4,) and (4,) into x and y coordinates.
We start of with 4 points, where we also have the z value (the function evaluation on the first two coordinates), so a shape of (4, 2+1)
coords.shape
We'll be receiving x and y values in the following shape, so this is what we have to work with:
_x = coords[:, [0]]
_y = coords[:, [1]]
_x.shape, _y.shape
Since the rotation_matrix is of shape (2,2), the final einsum operation is:
mi, ij -> jm
or
(nr_points, 2), (2, 2) -> (2, nr_points)
__x, __y = np.einsum('mi, ij -> jm', np.hstack((coords[:, [0]], coords[:, [1]])), rotation_matrix)
__x, __y
Again, let's wrap everything up in a single function to see where we are up until now.
plot_function_2d(himmelblau(), contour_log_scale=False)
Using both 2D and 3D plotting function we can combine them to be shown on the same figure, and using the same rotation.
As allways, we will first sketch the code in bulk then tie it up in a single function.
fig = plt.figure(figsize=(26, 10))
ax_3d = fig.add_subplot(1, 2, 1, projection='3d')
ax_2d = fig.add_subplot(1, 2, 2)
angle = 225
function = himmelblau()
plot_function_3d(function, ax=ax_3d, azimuth=30, angle=angle)
plot_function_2d(function, ax=ax_2d, angle=angle)
The 2D and 3D plots can now be combined in a unified function that you can see bellow:
from optimisations.functions import *
plot_function(himmelblau(), angle=225, contour_log_scale=False)
plot_function(mc_cormick(), angle=225, contour_log_scale=False)
plot_function(holder_table(), angle=225, contour_log_scale=False)
Let's see all of our functions defined
from optimisations.functions import Function
plot_all_functions(Function)